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Abstract - A particle in the attractive Coulomb field has an interesting property: its specific 
heat is constant and negative. We show, both analytically and numerically, that when a classical 
Hamiltonian system stays in weak contact with one such negative specific heat object, its statistics 
conforms to a fat-tailed power-law distribution with power index given by C/fcs — 1, where fcs is 
Boltzmann constant and C is the heat capacity. 



, Introduction. In 1968 Lynden-Bell and Wood |1 
"^pointed out an interesting fact: self-gravitating systems 
>have negative specific heats. As Lynden-Bell later on ex- 
plained [2] this fact seemed quite natural to astronomers, 
who know for example that when a star gains energy it ex- 
gpands and cools down 3|; but it appeared incongruous, if 
1 not completely wrong, to statistical mechanists, who learn 
' O from textbooks that specific heats are necessarily posi- 
gtive. The paradox was resolved by Thirring [4j: While it is 
true that any system in weak contact with a thermal bath, 
L — hence characterized by the canonical distribution, neces- 
sarily has a positive specific heat, isolated systems, charac- 
^ terized by the microcanonical distribution, may well have 
negative specific heats. Since then, many models showing 
CO microcanonical negative specific heats have been reported, 
"^see, e.g ., [5 -11 , their statistical mechanics has been dis- 
and negative specific heats have been ex- 



.cussed 12 13 



perimentally measured in thermally isolated small atomic 
clusters 



14p5] . Recently it has also been pointed out that 



, if a system is in strong coupling with its environment, like- 
• »wise it may display negative specific heat 16-23 . 
. ^ Inspired by Ref. [24] , here we consider an ordinary sys- 
^^tem S with Hamiltonian 7Js(x,p) and weakly couple it 
?— 1 via an interaction energy h(x, X), to a second system 
^with Hamiltonian Hc(X, P), possessing a constant neg- 
ative microcanonical heat capacity C < 0, i.e., 

ff(x, p, X, P) = £T s (x, p) + H C (X, P) + h(x, X) . (1) 

For the sake of clarity we recall that the heat capacity C is 
defined as the derivative of a system's energy with respect 
to the temperature, C = dE/dT, while the specific heat c 
is defined as the heat capacity per unit mass. 



The archetypical example of a system with constant 
negative microcanonical heat capacity is that of a single 
particle in the gravitational (or the attractive Coulomb) 
force field, for which C = —3/2 [2j. Throughout the paper, 
temperature is expressed in units of energy. In these units 
ks, Boltzmann's constant, is dimcnsionless and equal to 
1, and the heat capacity C = dE/dT is dimensionless as 
well. 

Our main result is that, provided the total system sam- 
ples the microcanonical ensemble, the system S samples 
the power-law distribution 



P(x,p) 



[H s (x,p) - E tot ] 



c-i 



J dxdp[i/ s (x, p)-E tl 



1C-1 : 



C < 0. 



(2) 



where E to t is the (conserved) energy of the total sys- 
tem. To express this result in Eq. ^ in the usual set 
of units, where the temperature is measured in Kelvin, 
and both the heat capacity and ks have the dimension of 
Joule/Kelvin, one should replace C in Eq. |2]) with the 
dimensionless ratio C /ks- 

We illustrate this result in Fig. [T] A neutral particle, 
with Hamiltonian Hs = p 2 /2m, is confined into a box and 
makes elastic collisions with a particle carrying a charge q 
subject to the Coulomb field generated by a fixed charge 
—q. The charged particle has constant negative heat ca- 
pacity C = —3/2, and its Hamiltonian reads: 



H, 



C=-3/2 



(X,P) = P 2 /2M -a/|X| 



(3) 



q 2 /Atteq with Eq the dielectric permittivity of 



where a 

vacuum. According to our main result, Eq. ([2]), since 
the charged particle has constant negative heat capacity 
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Fig. 1: (Color online) Illustrative example of our main results, 
Eq. (J2j . A neutral particle (orange) moves freely within a box, 
making elastic collisions with the box walls and a charged parti- 
cle (black) immersed in the attractive Coulomb field generated 
by a fixed point charge (red). Because the heat capacity of the 
charged (black) particle is C — —3/2, the velocity probability 
distribution function of the neutral (orange) particle conforms 
to a power law with power index —3/2 — 1 = —5/2, Eq. Q. 



C = —3/2, the velocity probability distribution function 
of the neutral particle is given by the following power law 



P(p) 



(p 2 /2m - 



£tot)- 5/2 



(4) 



Derivation of the result in Eq. ([2]). Assuming 
(i) the microcanonical distribution for the total system, 
and (ii) weak coupling h, the system S marginal distribu- 



tion p(x, p) can be written as 25} 27] 

n c [E tot -H s {x,p) 



P(x, p) 



where 



Q C (E) = J dXdP5[E - H C (X,P) 



(5) 



(6) 



is the density of states of the negative heat capacity system 
whose energy is E, and 

fitat(£tot)= /dXdPdxdp5[£^t-fl-(x,p,A:,P)] (7) 

is the density of states of the compound system. Here S 
denotes Dirac's delta function. 

We note that if a classical system has a constant negative 
microcanonical heat capacity C < 0, then its density of 
states is of the form: 



n c (E) cx (-E) 



c-i 



c < 0. 



(8) 



For the derivation of this density of states see Appendix 
A. In arriving at Eq. ^ we adopt the convention to set 
the zero of the energy of constant negative heat capacity 
systems as the lowest energy corresponding to unbound 
trajectories. Thus the density of states is only defined for 
negative E's, and it diverges for E > 0. 



Using Eq. Q in ^ one arrives at the result in Eq. 
([2]). It is important to stress that in Eq. ([2| the system 
S energy has a lower bound, which is conventionally set 
to zero, and that the total energy E tot must be negative, 
thus ensuring that the energy E — E to t — Hs of the system 
with negative heat capacity is negative at all times. This 
is necessary in order that the C-system always stays on 
bounded trajectories and never escapes to infinity. 

Numerics. — To corroborate our main result, Eq. pi, 
we simulated the dynamics of the system depicted in Fig. 
[TJ using a symplectic integrator [28]. We focussed on the 
probability density p(Eg) of finding the neutral particle 
(our S system) at energy Es for a total (negative) sim- 
ulation energy E to t- According to Eq. Q this is given 

by 

[Es — Etot 



P{E S ) 



-s/ 2 4 /2 

J °°dEs[E s -E tot ]- 5 / 2 E 

pl/2 



1/2 



(9) 



where the term E s derives from the density of states 

^s(Es) oc E s ^ 2 1 of the neutral particle. Note that for 
large E s , p{E s ) cx _E£ 2 . 

Following Ref. [27], we chose /i(x,X) = V LJ (\x-X.\), 
where Vlj is the truncated Leonard- Jones potential 



V LJ (r) 




+ e. 



\r\ > 2 1 / 6 a 
H < 2 x / 6 cr 



(10) 

We employ the same potential for the walls of the con- 
fining cubic box, in which the particle moves freely. In our 
simulation, we set e, a, and M (M is the mass of the par- 
ticle carrying the charge q) as the units of energy, length, 
and mass respectively. 

Following Ref. 27 , in order to avoid the singularity of 
the attractive Coulomb field at X = 0, in the simulation 



we employ the Plummer potential 29 : 

— n 

ip b (X) = 



Vx 2 



b 2 



(ii) 



Using this potential in place of the pure Coulombic po- 
tential <ys(X) = — a/|X|, induces a deviation of the den- 
sity of states from the pure power-law form f2_ 3 / 2 (£') = 
(-E)- 5 / 2 , resulting in a cut-off in the expected power law 
distribution p(Eg). This cut-off moves to higher energy 
values as b becomes smaller. 

Fig. [2] displays the result of our simulation for b — 
10~ 2 , Etot = " 1 5 a box of side length L = 10, a = 5, 
and m = y3- The computed energy distribution function 
for the neutral atom excellently agrees with the expected 
power-law form in Eq. ^ over three decades. At higher 
energies the distribution displays the expected cut-off due 
to the truncation in Eq. ( 11 ). Note that when the neutral 



atom has high energy the charged particle has low energy 
and stays close to the bottom of the Plummer potential, 
Eq. (11). The condition Etot < ensures that the energy 
E = Etot — Hs of the charged particle is also negative, 
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Fig. 2: (Color online) Energy probability density function 
p(Es) of a neutral particle performing Lennard- Jones colli- 
sions, Eq. (10) , with the walls of a confining box, and with 
a charged particle in the Plummer potential, Eq. (111. Black 



dots: numerical simulation. Red thick solid line: analytic cal- 
culation based on the Plummer potential, Eq. (11) . Black thin 
solid line: analytic calculation based on the pure gravitational 
potential <p(X) = Eq. |9|. Simulation parameters: 

E tot =-l,b= 10~ 2 , L = 10, a = 5, and m = y/3. 



thus ensuring that its motion remains confined at all times. 
With our choice of box side length L = 10, and E tot = — 1, 
the confinement keeps the charged particle within the box 
at all times. 

Methods. The simulation was performed by means of 

We 



an implicit Runge-Kutta symplectic integrator 28 



simulated the evolution of the whole system with a time 
step size At — 10 -5 and total number of time steps 
N = 1.28 x 10 12 . Note that the units of time are y/M Jea. 
In order to calculate the statistics, we divided the i?s-axis 
in 10 4 intervals of size AE$ = 0.0499, and counted how 
many times each interval was visited during the simula- 
tion. We adopted a sampling time r — 2 x 10 5 At. The 
resulting histogram was then normalized and divided by 
AEs to give an estimate of the probability density. Fi- 
nally, for better visualization in the log-log plot in Fig. 
[2j the histogram was coarse-grained with a series of in- 
tervals with exponentially increasing size, viz., we took 
AE% +1 /AE% = 1.1, where AEg is the size of the nth 
interval and AE S = AE S . 



Comparison with the Finite Bath Statistics. 

When C > 0, one finds the expression $lc(E) oc E C ~ Y 
for the densities of states, see Appendix A. This differs 
from Eq. Q by a minus sign in front of E. Note that 
for constant positive heat capacity systems, whose energy 
is bounded from below, we adopt the convention that the 
zero of the energy corresponds to this lower bound. Ac- 
cordingly, in this case one obtains for the probability den- 
sity function the same expression as in Eq. |2j), but with 



H s - E tot replaced by E tot - H s : 

[E^-Hsi^p)} - 1 



p(x,p) = 



/ dxdp[£ tot -ffs(x,p)] 



C-l ' 



C>0. (12) 



The statistics in Eq. (12) is sometimes referred to as the 
"finite bath statistics" 



20 



because it is induced by the 
weak coupling of the S system with a system with a finite 
positive heat capacity, like, for instance a finite collection 
of N free particles, for which C = 3N/2. As we will see 
below, this can also be achieved with one single particle in 
a properly chosen potential. The finite bath statistics in 
Eq. ( 12 ) interpolates between the microcanonical distri- 
bution p(x,p) oc 5[Eg — i?s(x, p)] and the canonical dis- 
tribution p(x, p) oc e~- Hs ( x,p )/ T . These two limiting cases 
are obtained when C approaches and +oo respectively 



26 30 31 . Note that when C > 0, E tot - H s in Eq. (12) 



is non-negative, as a consequence of the fact that the sys- 
tem S energy may not overcome the total (non-negative) 
energy -Etot- While the negative C distributions in Eq. 
([2]) have fat power-law tails, the positive C distributions 
in Eq. (12) exhibits sharp cut-offs at Hs — E to t- 



Implementing constant heat capacity systems 
with single particles. — In three dimensions, constant 
heat capacity systems can be obtained with a single par- 
ticle governed by the Hamiltonian: 



H c Qt, P) = P 2 /2M + gC\X/L\ 



6/(2C-3) 



(13) 



with constants g, L > with units of energy and length 
respectively, and C </ [0, 3/2]. As shown with explicit cal- 
culations in Appendix B, the associated density of states 
is proportional to E ^ 1 for C > 3/2, and is proportional 
to (— E) ^ 1 for C < 0. Accordingly, when a system 
with Hamiltonian Hs bounded from below stays in weak 
contact with one such constant heat capacity system, its 
statistics conforms with either Eq. ([2| or Eq. (12 1, de- 



pending on whether C<0orC>3/2, respectively. The 



potential in (13) includes many special cases of interest. 



For C approaching 3/2 the potential approaches a box 
potential, namely the system is a particle in a box, which 
knowingly has a constant heat capacity of 3/2. For C = 3, 
we have a harmonic oscillator. For C = —3/2 we have the 
1/r case of Eq. ([3]). For C — > ±00, the potential ap- 
proaches the logarithmic form: 



C|X| 6 A 2C - 3 ) ~31n|X| + <7, \C\ > 1 



(14) 



which, as we have shown recently 27 , produces the Gibbs 
distribution. Thus, not only does the canonical ensemble 
emerge from the coupling to a system with positive infinite 
heat capacity, but also from the coupling to a system with 
negative infinite heat capacity. The case C — is excluded 
from Eq. (13) because it corresponds to an unbound free 



particle, whose density of states, accordingly, diverges for 
all energies. For < C < 3/2, the exponent 6/(2C — 3) 
is negative, hence one must chose g < to ensure that 
the potential is attractive. The density of states diverges 
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nonetheless because it involves a divergent integral, see 
Appendix B. 



Similar expressions as in Eq. (13 1, but with a different 



exponent for |X| can be found in spaces of dimensions 
other than three. 

Conclusions. — We have shown, both analytically 
and numerically, that weak coupling to a system with 
constant negative heat capacity leads to the emergence 
of velocity probability distributions with fat power-law 
tails. The power-law emerges because the ^-system may 
withdraw indefinitely large amounts of energy from the in- 
finitely deep algebraic potential "well" that systems with 
constant negative heat capacity possess. 

Because of the fat tails it is possible that the negative 
C distributions may not be normalizable. This happens, 
for instance, when the S'-system density of states goes 
asymptotically as Eg, with a > — C. In this case the 
5-system "wins" over the C-system and the integral of 
^c[-Etot — Eg]Qg[Eg], which asymptotically goes like oc 
E°+ a ~ 1 , diverges. For a > -C + l, n c [E tot - E s }n s [E s } 
asymptotically increase with increasing energies, while for 
— C + 1 > a > —C it vanishes for Eg — > oc. As 
examples, consider the case of Fig. [T] but with many 
neutral particles. With two neutral particles we have 
-C + 1 = 5/2 > a = 2 > -C = 3/2, that is the dis- 
tribution would not be normalizable but would vanish for 
large Eg. With three or more neutral particles we have 
a > 7/2 > 5/2 = —C + l, and the distribution would 
be increasing for large Eg. Arguably the system never 
reaches equilibrium in such cases. 

For our example of a neutral particle making collisions 
with a particle in the 1/r potential field, Fig. [I] , the 
energy distribution has a power index —2, which is very 
close to the exponent observed in the energy distribution 
of cosmic rays 32p3[ ~ This suggests that the power index 
~ — 2 of the cosmic ray energy distribution might originate 
from multiple collisions of the cosmic particles with objects 
obeying 1/r potentials before they reach our instruments. 
We leave this as an open question. 
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Appendix A. Density of states for systems with 
constant heat capacity. — A classical system has a 
constant negative microcanonical heat capacity C if and 



1 Typically one looks at the cosmic ray differential spectrum 
32 33], which is the derivative of the probability distribution with 
respect to energy. The cosmic ray spectrum has a power index ~ — 3. 
Accordingly, the energy probability distribution has a power index 
~ -2. 



only if its density of states is of the form: 

n c (E) = fi-Ef- 1 , C<0 



(15) 



with some constant / > 0. 

To see this consider the microcanonical heat capacity: 



C(E) = [dT(E)/dE]~ 



(16) 



where T(E) is the microcanonical temperature. For con- 
stant C this implies 



T(E) = E/C + a, 



(17) 



with some integration constant a. Using the definition of 
microcanonical temperature 

T(E) = [dln^ c (E)/dE]- 1 , (18) 

one obtains, after a further integration 

$ C (E) = d{E/C + a] c , (19) 

with some integration constant d. Here 

§ C (E) = J dXdP9[E- H C (X,P)}. (20) 

is the phase space volume, and 9 denotes Heaviside's step 
function. Adopting the convention $0(0) = 00, sets a = 0, 
hence $ C (E) = d[E/C] c . Recalling that C < 0, this 
can be recast as $>c(E) oc (— E) c . Using the well-known 
relation 

n c (E) = d<£c(E)/dE, (21) 
we obtain Eq. (Il5| . 



To see that the reverse is also true, consider Eq. (21 1 



We have, using Eq. (15 I 

■ E 



®c(E) = [ n c {E')dE' = -f(-Ef/C . 

J —oc 



Using (18), we obtain: 



T{E) = E/C. 



(22) 



(23) 



From Eq. (16) it then follows that the heat capacity is C. 



Likewise, a classical system has a constant positive mi- 
crocanonical heat capacity C if and only if its density of 
states is of the form: 



n c (E) = fE 



C-l 



c > 



(24) 



with some constant / > 0. 

To see this the argument proceeds identically to the 
argument above, with the only difference that in the case 
C > the constant a is set to zero as a consequence of the 
convention §c(0) = 0. 

The argument supporting the reverse implication also 
proceeds identically as above, with the only difference that 
now the phase volume reads &c(E) = fE/C. 



p-4 



On the origin of power-laws in equilibrium 



Appendix B. Density of states for the Hamilto- 

We begin with the case C < 0, 



nians in Eq. (13) 



implying that the exponent 

7 = 6/(2(7- 



3) 



(25) 



lies in the interval — 2 < 7 < 0. Using Eq. (|20| with 



Eq. ( 13 1, recalling that E, C < 0, the phase volume reads, 



after integration over dP 
*ciE) = {2 ^^ J dX[E - gC\X\-<]^ 



(26) 



r(5/2) (-jO)»/i y 

where we switched to spherical coordinates in the sec- 
ond line and we performed the change of variables y = 
p(gC I E) 1 ' 1 in the third line. The limit of integration 
(E/gC) 1 ^ 1 represents the turning point of the closed or- 
bit of negative energy E. In the interval —2 < 7 < the 
integral converges and can be expressed in terms of the 
gamma function T(x), as: 



dyy 2 [y 1 - 1] 



3/2 



^ T(-C) 
4 r(3/2-C) 



(27) 



Note that 3/2 + 3/7 = C. Hence $ c oc {-E) c . Accord- 
ingly, using Eq. (21 ), we find that the density of states is 
of the form tt c cx (-E) ' 1 . 

For 0< C < 3/2 it is 7 < —2, hence one must now chose 
g < in order that the potential is attractive. One then 
arrives at the same expression as for the case C < 0. In 
this case however the integral J Q dyy 2 ^ 1 — l] 3 / 2 diverges. 

For C > 3/2 it is 7 > 0. One can proceed as above, 
recalling that now E, C > 0, to obtain 



(2irM) 3 / 2 r^/aC) 1 '- 1 
{ ] 1 dpp\E-gCp^ 



47T 



r(5/2) Jo 

(2^M)^J 1 dyy 2 [l-y^/ 2 c 
r(5/2) {gCfh 



(28) 



The integral in the numerator is well behaved for all 7 > 0, 
and can be expressed as: 



Jo 



fRTjG -1/2) 
4 r(C + l) 



(29) 



Using Eq. (21 ), we find that the density of states is of the 
form Q c ocW 1-1 . 
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